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Abstract 

This study of occupation time densities for continuous-time Markov processes was inspired 
by the work of E.Nir et al (see [13]) in the field of Single Molecule FRET spectroscopy. There, a 
single molecule fluctuates between two or more states, and the experimental observable depends 
on the state's occupation time distribution. To mathematically describe the observable there 
was a need to calculate a single state occupation time distribution. 

In this paper, we consider a Markov process with countably many states. In order to find a 
one-stete occupation time density, we use a combination of Fourier and Laplace transforms in 
the way that allows for inversion of the Fourier transform. We derive an explicit expression for 
an occupation time density in the case of a simple continuous time random walk on TL. Also we 
examine the spectral measures in Karlin-McGregor diagonalization in an attempt to represent 
occupation time densities via modified Bessel functions. 

Primary Subjects: 60J35, 47N30, 47N20, 33C90, 33C10 

Keywords: occupation time, Bessel functions, FRET spectroscopy, orthogonal polyno- 
mials 

1 Introduction 

The occupation time densities for continuous-time Markov processes that live on countable 
state space was a subject of intense research in the 60s, 70s and early 80s. We would like 
to refer the reader to [5], [6], [H], [H], [16] and [3j (many published in the Journal of 
Applied Probability) for some of the results in the field. With the exception of [3], the 
main instrument was the multidimensional Laplace transform. Also, the reader can check 
[T3] and [2] for the most recent developments in the field. 

This paper is a mathematical follow-up to the research done by E.Nir et al [13] in the 
field of Single Molecule Fluorescence Resonance Energy Transfer (FRET) spectroscopy, 
where a single molecule fluctuates between two or more states, and the experimental 
observable depends on the state's occupation time distribution. While working on [13] 
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the authors have noticed that the single state occupation time densities, when computed 
via randomization technique (i.e. multiple infinite sums) can often be represented via 
modified Bessel functions of the kind 



E 



1 / z\2k+n 



V 



k\{k + n)\ \2) 

In this paper, we use spectral theory in an attempt to find a rigorous explanation for the 
relationship between occupation times and Bessel functions. We will show a connection 
between a spectral measure of a generator and a Laplace transform of a single state 
occupation time distribution taken with respect to a time variable t. 

The occupation times for birth-and-death chains were studied by Karlin and McGregor 
with orthogonal polynomials in [TO] following the paper of Darling and Kac, see [1] . Both 
papers considered occupation times for Markov processes when t is taken to cxd, while this 
paper concerns with explicit expressions for a given fixed time interval [0, t]. 



2 Approach and results 

In this section we will state and prove theorems and formulas that will relate occupation 
times to the spectral measure of a generator, we will find an explicit expression for a 
one-dimensional symmetric random walk (see Theorem 2.2), and find a representation of 
occupation time distributions via modified Bessel functions. 



2.1 General case: spectral representation 

Consider an irreducible discrete Markov process with the generator matrix (or operator, 
if the state space is infinite). 



Q 



•^1,0 



Ao.i 



...\ 
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•••/ 



Let {0, 1, ... } be a countable state space, and for a time interval [0, t], let /^(t, x) denote 
the probability density function for the occupation time associated with state 0, given 

1 



that the continuous-time process commences at state k. Denote cq 








First, we derive the following expression for the occupation time density fo{t,x). 
Theorem 2.1. The Laplace transform w.r.t. time t of fo{t,x) can be written as 

1 I X 



Lf,{si,x) 



Sih{s^ 



exp 



'h{s^] 
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where h{s) — — {{Q — si) ^Cq, eo). 

Proof. The integral equations relating {/fe(t, x)}fc=o,i,... can be produced via conditioning 
as follows: 

foit.x) = e-(^- -^i^°''")*5t(a;) + ^ j fk{t - y,x - y)Xo,ke~'^^"^- "^^^^^'-^''^dy, 
fj{t,x) = e~^^^- ^^i^^'"''>^5Q{x) + / fk{t - y,x)Xj^ke~^^"'- ""^'^''""^^dy 



k: k^j 

forj = l,2,... . 



Plugging in ip — t — y into the above equation, we obtain 



t 

-^J,m)(t-l/') 



k: k^j 

for J = 1,2,... . 
Now taking the Fourier transform with respect to x we arrive at 

COD ft 



dip 



/OO /"f 
/ fk{'ip,x)Xj,ke~^^^- "^^^^^'-^^^^-^Uipe^'^^'dx 
„. -oo Jg 



k: k^j ' 

forj = l,2,... 



The above equations simplify to 

k: fc^C^^ 



k: k^j 

for j = l,2,... . 



Differentiating w.r.t. variable obtain 



( Ao,m -«S2 j /o(^,S2) + ^/o(^,S2) = Ao,jk/fc(^, S2) , 

( ^^'"^ ) ^2) + ^2) = XI ^hkh{t, S2) (i = 1, 2, . . . ) 

\m: m^j / k: k^j 
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Observe: /^(O, S2) = 1 for all j. Our next step is to take the Laplace transform w.r.t. 
variable t: 



1,^2 



si,S2) = 1+ ^ Aj-fcL^;(si,S2) (i = l,2,...)- 



The above system of equations can be rewritten via the spectral decomposition of the 

J > -, . So we 

proved the following spectral identity 

" L^,(Si,S2) 









" 1 " 


enerator operator Q as follows. Let Lj{si,S2) = 


L/,(S1,S2) 


and 1 = 


1 



{Q - siI)Lj:{si, S2) = -1 - is2 








Thus L^-(si, S2) = -(Q - si/) ^1 - is2Lj^^{si, S2){Q - sil) ^Cq and 

Lf^{si,S2) = - ((Q - si/)~4,eo) -is2Lj:^^{si,S2) {{Q - si/)"^eo,eo) • 
Therefore the Laplace-Fourier transform of /o can be represented as 

-((Q-si/)-il,eo) 



L;„(si,S2) 



I + is2 HQ - sil) ^eo,eo) 



(3) 



Observe that Lj:^(si,0) = /[o+oo)/ir^ '^^^ fk(t,x)dxdt = ^ for all k. Substituting S2 = 
into gives 

^-{Q-siI)l = -l 



Sl 



which is obviously true. This also implies 



iQ-siI)-'l = 1 . 



Therefore ([3]) can be simplified to 

L/o(Sl,S2) 



1 - is2h{si) 
4 



(4) 



1 . 1 VU V UliCgU V , IN .iVlClCUlLli aiiU J-L/.iNll 



vyuuupaiiuii iiiiico via jjcsaci iuiiuiiuiia 



where h{s) = - ((Q - sJ)-^eo, Cq) = {{J^ e-'^e^^dt)eo,eo) = J^^ e-'%{0,0)dt. The 
Fourier transform can be inverted via complex integration over a lower semi-circle contour 
with the radius converging to infinity: 

L/o(si,a;) = — r^expl— — = — — ]— 7 exp I ^— -I. (5) 

sih{si) [ h{si) ) si{{Q - sil) ^eo,eo) [{{Q - sil) ^eo,eo)J 

□ 

2.1.1 Example: Two-state Markov processes 

Consider a two-state Markov process with generator Q = ( I • Then 

V J 

(D- r^i-i - ~^ ( Ai + '^l ^ \ 

sl + {\ + ^,)s^\ /X \ + s,) 

and (P) implies 

Si + /i 

Now, formula (29.3.81) of [1] gives us the following Laplace transforms 

/ /o(2v/at)e-P*dt = -et and / —lA2^t)e-P^dt = —{ep - 1), 
Jo P Jo Vt Va 

where Iq{-) and /i(-) are modified Bessel functions. Next, we rewrite the above identities 
as follows 

.1 a 



P 

and 



OO 



/ h{2^a{t-x))e-P'dt. 
\/t - X 



Let a = A/ix and p = Si + /i. Plugging in, we get 

1 f > \ ^^J: 1 a 



Sl+ n 



-(«i+M)Xe^ = e-P"-ep = / Io{2^Xfix{t - x))e-^''e-'''dt 

P J X 



and thereforethe inverse Laplace transform of -^^—e'^^'^^^^'c^+p- is 



^e-Axg-^(t-x)j^^2^A/ia;(t - x)) 



for < X < t. 
Similarly 



e 
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which can be rewritten as 



e 



Jx y/t-X Jo 



Xfix 



Therefore, the inverse Laplace transform of e ^(*i+'^)en+f' is 



Here we do not divide by zero when x = t as the \/t — x cancels on top and the bottom. 
Adding the terms together, we obtain 



^-\x^-ii{t-x) 



fo{t,x) = e-^*5t(x)+Ae-^^e-'^(*-^)/o(2v/A/ix(t - x)) + ^^h{2^ Xf,x{t - x)) 
for < X <t. 

The above equation was originally derived in [T3j via two-dimensional Laplace trans- 
form. One can also derive it via randomization, where the infinite sums are easily recog- 
nized to be the corresponding modified Bessel functions. 

2.1.2 Example: Three state Markov chain. 

Here, the inversion of Lf^{si,x) as expressed in ([T]) with 

1 7i 72 



h{z) z + pi Z + /32 

can be expressed via convolutions of modified Bessel functions. 

The last two examples prompted us to look closely at the structure of the occupation 
time densities, and in an attempt to understand the mechanics of decomposing them via 
cylindrical functions /n(-)- 

2.2 Continuous time birth-and-death chains and related pro- 
cesses. 

In the case of a birth-and-death process, the spectral representation ([1]) of Ljg(si,x) can 
be expressed via orthogonal polynomials. 



Let ttq = 1 and TCk = ° ^ for all k > 1. Observe that vr = [vrfO), 7r(l), 7r(2), . . . 



AoAi...Afc_ 

satisfies the detailed balance (reversibility) condition for the process: 

TTfcAfc = TTfe+i/ife+i for all A; > 0. 

Let Po{s) = 1, -Pi(s), P2(s),... (where each Pk{s) is a polynomial of kth. degree) be 
constructed recursively as the coordinates of an eigenvector 



P[s] 



P2{S) 



satisfying {Q — sI)P[s] = . 



6 



1 . 1 VU V UliCgU V , IN .iVlClCUlLli aiiU J-L/.iNll 



vyuuupaiiuii iiiiico via jjcsaci iuiiuiiuiia 



We recall the results of [H] and [H], where it was shown (extending a theorem of J.Favard) 
that there is a probability measure n on (— cxd, 0] with infinite support such that the 
polynomials {Pk{s)}k=oA,... are orthogonal w.r.t. measure /x, 



Pkis)Pmis)dlJ,{s) 
(-00,0] ^fc 



(6) 



The expression ([T]) for Lj-u(si,x) can be rewritten via the Cauchy transforms w.r.t. 
the spectral probability measure /i as follows. 

{{Q-sJ)-'ek,eo)= [ ^^dfi{x) = CiPkdfi){si), 

J{-oofl] X — Si 

where C{gd^){s) = /(_^ qj ^^dfi{x) denotes the Cauchy transform (w.r.t. dfi) of g. Here 
the function h{s) of can be expressed as 

h{s) = -[ 

2.2.1 Birth-and-death process with equal rates 

Here, we will compute the occupation time density for a birth-and-death process with 
forward rates A = 1 and reverse rates = 1, i.e. the process whose generator is a simple 
Jacobi operator 

\ 
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where r > 0. 



Theorem 2.2. The zero-state occupation time density for a hirth- and- death process with 
equal rates can he expressed via modified Bessel functions as follows: 

/o(t, x) = e-''6oit -x)+ re^^-'^'^'h (2^J [t - x){t + {r - l)x)) ■ 

+ _e(^-)-^Vi (2^it-x)it+ir-l)x)) ■ 1|.<,| . 



Proof. Equation ([2]) translates as 

L^(Si,S2) = — — - — — + — — - — — L^(Si,S2) 
■'° r + Si — r + Si — iS2 ■' 

H(-'uS,) = ^ + ^L;^_.(-l'^2) + ^L;^,,(.l,32) (A: = 1, 2, . . . ) 



where L^^(si, S2) again denotes the Laplace transform of Fourier transform of fk. In this 
recurrence relation, let 



Si 



Then Ik satisfy the following recurrence relation, 



/fc_i(si,S2) + TTT — h+iisi,S2) (A; = 1,2, . . .), 



2 + si ^' 2 + si 
Solving the characteristic equation, 

- (2 + + 1 = 

get 



k{si,S2) = Qil(Si,S2) 



2 + si + V^i + 4si 



+ Q!2(Sl, S2) 



2 + Si - ^/4T4^ 



Observe that L^^(si,0) = /[o,+oo)/r^ ^^^fk{t,x)dxdt = ^ and 



L/^(si,S2)= / / e 

'[0,+oo) 

That is 

Hence, since si > 0, 



x)dxdt 

'[0,+oo) 

-S2) — as /c — > 00 



-S\t+iS2X 



So{x)dxdt — — as A; ^ 00 . 



Sl 



lk{si,S2) = k{si,S2) 



2 + si- ^sfTI^ 



Now the first recurrence relation for can be rewritten as 



Therefore 



— + — 

r + SI-IS2 r + Sl - IS2 



(L;„(.i,.2)--) 



1 \ 2 + Sl - + 4si 1 



+ 



(2-r)si + rv/s? + 4si 



2si s2 + |((2-r)si + rVs? + 4si) ' 
Once again, using complex integration, we arrive to 

{2-r)sl + r^/sl + Ail 



L/o(si,x) 



2si 



exp<{ -- ( (2-r)si + r^/s^ + 4si 
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Observe that one can use the same characteristic equation in order to find the expression 
for the corresponding orthogonal polynomials: 



k 



where the spectral measure will satisfy 



(-00,0] 



X 



s (2 - r)s + r\/ + 4s 



We will now invert the Laplace transform by decomposing Lj(,(si,x) as follows 

L/o(si, = + \Vn + 2rVni , 



where 



Vi = exp {-|(2 - r)si| ■ exp |-|r-^s? + 4si| , 
Vn = ;^7^1=^^P {-|(2 - • exp {-|rv^s? + 4si} 



and 



Vui = ^7=^ {-f (2 - r)s^} ■ exp |-|r^s? + 4s,} . 
We will quote a Laplace transform formula (29.3.91) in [T|: 

/ e-^*e-2«*/Q -aVt^ - dt = {k > 0) . 

Jk \2 J A/s(s + a) 

First we will find the inverse-Laplace transform of Vm- Taking s = si, a = 4 and k 
in (29.3.91) of ^, we get 



rx 
2 



e-^^V2*/nl2./t2_(l^) \dt 



^2/ / Vsi(si + 4) 



Multiplying both sides of the above equation by exp |_ (^ J^'''^ si|, and changing the vari- 



able to t := t + obtain 



2 



ViH = £ e-^^*e(2"'^)^-2*/o (2v/(t-a;)(t + (r-l)a;)) dt . (7) 
Therefore, the inverse of Vm is 

/:-'(7'm) = e(2-^)^-2*/o (2V(t-x)(t + (r-l)x)) ■ l{.<t}. 
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We differentiate -k? and integrate by parts in ([7]): 



Hence 



d_ 
di 



Vii - e-™e 



rx — Six 



C~\Vn) = e-'^'Soit - x) - 2e(2-'-)-2*Jo (2 ^(t - + (r - l)x)) ■ l{,.<t| 

+ ^=£=^^^==e(-^)^--/, (2v/(t-x)(t + (r-l)x)) ■ 1,.,,,. 
A/(t — + (r — ^ ^ 

In order for us to invert Vj, we will need (29.3.96) of [T], that states the following 

e-^^-^L=h (aVW^) dt = e-'^vW _ ^-ks^ > 0) . 
Here we let s = si + 2, a = 2 and k = thus obtaining 



Once again changing the variable to t := t + get 



foo {2-r)x-2t 



+ (r - l)x) 



h {2^{t-x){t+ (r- rft 



and 

Ii-\V,) = e-^%it -x)+ ==h (2v/(t-x)(t + (r-l)x)) ■ 1|.<,} . 

^y (t X jit -\- yv LJX J ^ 

We add up all three terms together, thus proving the theorem. □ 



2.2.2 How the occupation time densities are expressed via modified Bessel 
functions, !„.{■) 1 and the moments of the spectral measure 

Let mo, mi, . . . denote the moments of the spectral measure n, i.e. 

mj = / {—xy dfi{x) . 
J {-00,0] 

Consider a case where the spectral measure fi has bounded support, say supp{fi) C 
[—K, 0]. Then, for any z G C \ (—00, 0] such that \z\ > 2K, 



1 



h{z) 1 — miz~^ + m2-2^^ — . . . 
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1 + ^{miz^^ - m2Z~^ + ...y 



k=i 
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= z + nil — (^2 "~ Tn\)z ^ + 0(2; ^ . 
Recall ([T]). Now, we will consider the inverse Laplace transform of a function Fo(t, x). 
whose Laplace transform IjFq{si, x) = exp | — /^(f^y |- Here for an a > and < x <t, 



Fo(t, x) = — / - exp <( zt - } dz 



2''r^ J a—ioo ^ 



= - — : / - exp \^z{t — x) — rriix + {m2 — 'rnl)xz^^ — X(l){z^^)z^'^^ dz . 

J a-ioo ^ 

Here e-^<^(^~')^~' = 1 + ^^=2 ^^(2^)^"^- Thus 



Fo{t,x) = 6-^"^% (2y/(m2-m2)(t-a;)x 



+ e" 



fc=2 



t — X 



(m2 — m^)x 



4 2W (ma - mf)(t - x)a; 



by (29.3.81) in [T]. 

Another derivation of the same formula may allow further simplification. Since Q is 
the generator for a reversible Markov process, its spectrum is entirely contained inside 
(—00,0]. Here the spectrum is a subset of [— i^', 0]. Let m{z,x) = 2^^~^^^^ 
for s G (— cxD,0], v{s,x) = m_(s,x) — m+(s, x), where m_(s,x) = limsiom{s — ie,x) 
and m+{s,x) = limej^o^('S + ie,x). Observe that 'm{z,x) as \z\ 00. Solving the 
Riemann-Hilbert problem via Plemelj formula, obtain 



m(z, x) 



27Ti 



v{s, x) 



(-00,0] 



z — s 



ds for zeC\[-K,0]. 



Now, 



Fo{t,x) 



a+ioo 



miz, X 



exp \^z{t — x) — rriix + (m2 — m\)xz dz 



(-00,0] 



1 

2tH 



a—ioo 
a+ioo 



Z — S 



exp \^z{t — x) — rriiX + (m2 — m\)xz dz 



f (s, x)ds 



-mix 



(-00,0] ^^Q 



a+too 



2vr'j .1 a-ioo ^ 



k+1 



exp \^z{t — x) + (m2 — m\)xz ^} dz 



v{s, x)ds 



'^""Y^ / s^v(s,x)ds , , „, 
^A-oo,o] \{m2-mf)x 



t — X 



Iki'^y {irL2 - mf){t - x)x 



by (29.3.81) in [Ij. 

Observe that since (j){z)z~'^ has no poles at 00, for k > 1, 



s^v{s,x)ds=i z^-^e-''"^^'''^'~'dz 

(— cxD,0] ^7+U7_ 
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where 7+ connects the origin to 00 right above the negative half-hne, and 7_ connects 
cxo to the origin barely bellow the negative half-line. Observe that for k = 1 the above 
integral is zero. 



Observe that in general, an argument in Deift [7j (Section 2) shows that if a process is 
time reversible (i.e. satisfies a detailed balance condition) with a bounded generator then 
there exist a unique (spectral) probability measure ^ with compact support supp{fi) C M 
such that 



As a conclusion, let us list some of the open problems and directions for further research 
the authors are working on. 

• What properties of dfi would allow for the inversion of the Laplace transform Ljq(si, x) 
to be expressed explicitly via modified Bessel functions /„? 

• Exploring occupation time densities for a wider class of time reversible stochastic 
processes. 

• Interpreting reinforced processes studied in Kovchegov [12j as occupation time driven 
processes. 

Acknowledgment 

The authors wish to thank R.Burton and M.Ossiander for sharing thoughts on the subject 
of this paper. 

References 

[1] M.Abramowitz and I. A. Stegun (Eds.), Handbook of mathematical functions with 
formulas, graphs, and mathematical tables U.S. Department of Commerce (1972) 

[2] M.Bladt, B.Meini, M.F.Neuts and B.Sericola, Distributions of reward functions on 
continuous time Markov chains in Matrix-analytic methods: theory and application 
(Eds. G. Latouche and P. Taylor) (2002), pp.39-62 

[3] L.Bondesson, On occupation times for quasi-Markov processes J.Appl.Prob., 18, 
(1981), pp.297-301 

[4] D. A. Darling and M.Kac, On occupation times for Markoff processes Transactions of 
AMS, 84, (1957), pp.444-458 

[5] J.N.Darroch, Identities for passage times with applications to recurrent events and 
homogeneous differential functions J.Appl.Prob., 3, (1966), pp. 435-444 



Here fo{t,x) 




o(t,a;). 



3 Conclusions 




12 



vyuuupatiuii tiiiicB via ucssci luiiuLiuiis 



[6] J.N.Darroch and K.Morris, Passage-time generating functions for continuous-time 
finite Markov chains J.Appl.Prob., 5, (1968), pp. 414-426 

[7] P.Deift, Orthogonal Polynomials and Random Matrices: A Riemann-Hilbert Ap- 
proach Amer. Math. Soc, Providance, RI, (2000) 

[8] S.Karlin and J. L. McGregor, The differential equations of birth and death processes 
and the Stieltjes problem Transactions of AMS, 85, (1957), pp. 489-546 

[9] S.Karlin and J. L. McGregor, The classification of birth and death processes Transac- 
tions of AMS, 86, (1957), pp.366-400 

[10] S.Karlin and J. L. McGregor, Occupation time laws for birth and death processes Proc. 
4tli Berkeley Symp. Math. Statist. Prob., 2, (1962), pp.249-272 

[11] L.M.Kovaleva, On the occupation time in a given state for the simples semi-Markov 
system Teor. Veroyat. i Mat. Stat., 1, (1970), pp.100-107 

[12] Y.Kovchegov, Multi-particle processes with reinforcements Journal of Theoretical 
Probabihty, 21, (2008), pp.437-448 

[13] E.Nir, X.Michalet, K.Hamadani, T.A.Laurence, D.Neuhauser, Y.Kovchegov and 
S.Weiss, Shot-noise limited single-molecule FRET histogram: comparison between 
theory and experiments Journal of Physical Chemistry B, Vol.110, No. 44 (2006), 
pp.22103-22124 

[14] P.J.Pedler, Occupation times for two-state Markov chains J.Appl.Prob., 8, (1971), 
pp.381-390 

[15] B.Sericola, Occupation times in Markov processes Stochastic Models, 16, (2000), 
pp.479-510 

[16] W.S.Hsia, The joint probability density function of the occupation time of a three-state 
problem J.Appl.Prob., 13, (1971), pp.57-64 



13 



